Program main

!      subroutine Evaluate_Capitation_Payment (InsurSize, InsurClaimSize, InsurStdDev, InsrCapPayment, ProvrSize, ProvrProfitGoal, ProvrStdErr, ProvRiskAdjCapCharge)   

!****************************************************************************************
!
!! Evaluate_Capitation_Payment evaluates the capitation payment offer of an insurer
! in terms of whether it is actuarially adequate for a Health Care Provider.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    21 May 2013
!
!  Author:
!
!    Original FORTRAN90 version by Thomas Cox
!
      implicit none

	real InsurSize, InsurStdDev, InsurClaimSize, InsrCapPayment
	real ProvrSize, ProvrStdDev
 	real ProvRiskAdjCapCharge, ProvrCost, ProvrProfitGoal
      integer*4 i, j, k , l, m, n 
      character * ( 100 ) phrase
!
! Get information on insurer's claims
!

write (*, '(a)') 'Enter the size of the insurer, the insurers average claim size, the insurers claim size standard error and the insurers capitation payment offer:'
read (*, 4f12.0) InsurSize, InsurClaimSize, InsurStdDev, InsrCapPayment 

!
! Get information on Providers Operations
!

write (*, '(a)') 'Enter the size of the provider, the providers average claim cost and the providers profit goal expressed as an integer % of capitation revenues'
read (*, 3f12.0) ProvrSize, ProvrCost, ProvrProfitGoal

!
!	Calculate the provider's standard error, the provider's standard error
!	and assume the provider wants to meet/exceed their profit goal with 
!	probability = 0.9500
!
!	The provider's capitation payment is the solution to the equation:
!
!!	Prob[InsrCapPayment .GE. (1.0 + ProvrProfitGoal/100.0) * ProvrCost]
!
!	Where this probability is the value of 1 - CDF for the normal distribution at
!	the point 1.96 provider standard errors above the average claim cost 
!
!	Providers fail to meet their profit goals when the actual aggregate claims amount for 
!	their roster of patients exceeds their aggregate capitation payments and their 
!	profit goal
!
!
!	To earn their profit goals a providers costs must be less than 
!	(1.0000 - ProvrProfitGoal/100.0)
!
!	If this is not true the provider cannot attain their profit goal
!

	If (ProvrCost .GT. (1.0000 - ProvrProfitGoal/100.0) * InsrCapPayment) then
		write(*.'(a)") 'Costs exceed the level likely to yield desired profits. You may earn such profits in very good years, but not in average years'

	ProvrStandardError = InsurStdDev / (ProvrSize ** 0.5)

	ProvrCost + 

	ProvRiskAdjCapCharge = 
!
!	
!
!
!	ProvRiskAdjCapCharge = (1 + ProvrCost + , ProvrProfitGoalP
!
!	

write (*, '(a)') 'Profit/(Loss) %     Probability'
write (*, 2f12.4)') ProvrProfitLevel, ProvrProfitProbability
do i = 1,20



      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'Claims Generator For Risk Transferring Health Care Finance Mechanisms'
      write ( *, '(a)' ) '  FORTRAN90 version'


      CALL init_random_seed()         ! see example of RANDOM_SEED

      CALL RANDOM_NUMBER(r2)

OPEN(UNIT=12, FILE="GammaRNsout.txt", ACTION="write", STATUS="replace")
OPEN(UNIT=13, FILE="GammaRNsout2.txt", ACTION="write", STATUS="replace")


stop
end

subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  character ( len = 8 ) date
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  character ( len = 10 ) time
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end


      subroutine normal_01_cdf_values ( n_data, x, fx )

c*********************************************************************72
c
cc NORMAL_01_CDF_VALUES returns some values of the Normal 01 CDF.
c
c  Discussion:
c
c    In Mathematica, the function can be evaluated by:
c
c      Needs["Statistics`ContinuousDistributions`"]
c      dist = NormalDistribution [ 0, 1 ]
c      CDF [ dist, x ]
c
c  Modified:
c
c    24 March 2007
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c  Parameters:
c
c    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
c    first call.  On each call, the routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    Output, double precision X, the argument of the function.
c
c    Output, double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 17 )

      double precision fx
      double precision fx_vec(n_max) 
      integer n_data
      double precision x
      double precision x_vec(n_max) 

      save fx_vec
      save x_vec

      data fx_vec /
     &  0.5000000000000000D+00, 
     &  0.5398278372770290D+00, 
     &  0.5792597094391030D+00, 
     &  0.6179114221889526D+00, 
     &  0.6554217416103242D+00, 
     &  0.6914624612740131D+00, 
     &  0.7257468822499270D+00, 
     &  0.7580363477769270D+00, 
     &  0.7881446014166033D+00, 
     &  0.8159398746532405D+00, 
     &  0.8413447460685429D+00, 
     &  0.9331927987311419D+00, 
     &  0.9772498680518208D+00, 
     &  0.9937903346742239D+00, 
     &  0.9986501019683699D+00, 
     &  0.9997673709209645D+00, 
     &  0.9999683287581669D+00 /
      data x_vec /
     &  0.0000000000000000D+00,   
     &  0.1000000000000000D+00, 
     &  0.2000000000000000D+00, 
     &  0.3000000000000000D+00, 
     &  0.4000000000000000D+00, 
     &  0.5000000000000000D+00, 
     &  0.6000000000000000D+00, 
     &  0.7000000000000000D+00, 
     &  0.8000000000000000D+00, 
     &  0.9000000000000000D+00, 
     &  0.1000000000000000D+01, 
     &  0.1500000000000000D+01, 
     &  0.2000000000000000D+01, 
     &  0.2500000000000000D+01, 
     &  0.3000000000000000D+01, 
     &  0.3500000000000000D+01, 
     &  0.4000000000000000D+01 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        x = 0.0D+00
        fx = 0.0D+00
      else
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end


